#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue May  2 18:15:26 2023

@author: jcfq2
"""

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np
import shapely.geometry as sgeom
from cartopy.feature import ShapelyFeature


# import cartopy
import cartopy.crs as ccrs

from datetime import datetime

today = datetime.now()

def get_values_now():
    direct = today.strftime("%d%b%y").lower()
    time = today.hour*60*60+today.minute*60+today.second
    jd,jc,js,jnpd,jnpa = get_values(direct,time)
       
    return jd,jc,js,jnpd,jnpa


def get_values(direct,time):

    # time = seconds after 00:00 that day

    jupiter_rotationrate=360./(9.9258*60*60) #degree/secs (synodic)
    date = []
    diam = []
    cml = []
    sel = []
    npd = []
    npa = []
    
    f = open('jupiter_values.txt')
    for line in f.readlines():
        date.append(line[10:12]+str.lower(line[6:9])+line[3:5])
        diam.append(line[23:29])
        cml.append(line[33:38])
        sel.append(line[44:51])
        npd.append(line[55:65])
        npa.append(line[65:74])
    f.close()
    
    # date=np.array(date)
    diam=np.array(diam).astype(float)
    cml=np.array(cml).astype(float)
    sel=np.array(sel).astype(float)
    npd=np.array(npd).astype(float)
    npa=np.array(npa).astype(float)
    
    delta_cml=time*jupiter_rotationrate
    
    str_loc=int( np.where( np.char.find(date,direct) > -1) [0][0] )

    return diam[str_loc],(cml[str_loc]+delta_cml)%360,sel[str_loc],npd[str_loc],npa[str_loc]


# In[4]

class jup_night(ShapelyFeature):
    def __init__(self,              latit,longit,delta=0.1,refraction=0.01,
                   color="k",alpha=0.5,**kwargs):
        pole_lon = longit
        pole_lat = latit
        if latit > 0:
            pole_lat = -90 + latit
            central_lon = 180
        else:
            pole_lat = 90 + latit
            central_lon = 0
    
        rotated_pole = ccrs.RotatedPole(pole_latitude=pole_lat,
                                        pole_longitude=pole_lon,
                                        central_rotated_longitude=central_lon)
        npts = int(180 / delta)
        x = np.empty(npts * 2)
        y = np.empty(npts * 2)
    
        # Solve the equation for sunrise/sunset:
        # https://en.wikipedia.org/wiki/Sunrise_equation#Generalized_equation
        # NOTE: In the generalized equation on Wikipedia,
        #       delta == 0. in the rotated pole coordinate system.
        #       Therefore,              the max/min latitude is +/- (90+refraction)
    
        # Fill latitudes up and then down
        y[:npts] = np.linspace(-(90 + refraction),              90 + refraction,              npts)
        y[npts:] = y[:npts][::-1]
    
        # Solve the generalized equation for omega0,              which is the
        # angle of sunrise/sunset from solar noon
        # We need to clip the input to arccos to [-1,              1] due to floating
        # point precision and arccos creating nans for values outside
        # of the domain
        arccos_tmp = np.clip(np.sin(np.deg2rad(refraction)) /
                             np.cos(np.deg2rad(y)),              -1,              1)
        omega0 = np.rad2deg(np.arccos(arccos_tmp))
    
        # Fill the longitude values from the offset for midnight.
        # This needs to be a closed loop to fill the polygon.
        # Negative longitudes
        x[:npts] = -(180 - omega0[:npts])
        # Positive longitudes
        x[npts:] = 180 - omega0[npts:]
    
        kwargs.setdefault('facecolor',              color)
        kwargs.setdefault('alpha',              alpha)
    
        geom = sgeom.Polygon(np.column_stack((x,              y))) #polygon object
        # print(geom.shape)
        return super().__init__(
            [geom],              rotated_pole,              **kwargs)



jd,jc,js,jnpd,jnpa = get_values_now()
# jd,jc,js,jnpd,jnpa = get_values('11sep98',9*60*60)

print(jc)

grod_nio=np.array([80.23,81.36,82.14,82.69,83.04,83.17,83.09,82.80,82.13,80.24,77.25,70.63,65.35,60.44,55.22,51.35,49.84,50.40,52.30,54.68,57.01,59.32,61.49,63.72,66.07,68.08,69.48,70.45,71.12,71.70,72.29,73.09,74.17,75.58,77.24,78.85,80.23])
grod_sio=np.array([-62.72,-62.20,-61.79,-61.24,-60.45,-59.56,-58.96,-58.94,-59.69,-61.31,-63.76,-66.70,-69.82,-72.73,-75.16,-77.05,-78.45,-79.39,-80.05,-80.43,-80.60,-80.55,-80.30,-79.89,-79.11,-78.18,-76.92,-75.42,-73.71,-71.91,-70.13,-68.43,-66.90,-65.56,-64.40,-63.45,-62.72])

grod_n30=np.array([88.04,88.22,88.33,88.40,88.41,88.40,88.33,88.21,88.05,87.81,87.42,86.85,85.97,84.55,82.13,63.00,57.43,56.06,56.88,58.94,61.37,64.10,67.05,70.45,73.76,76.58,78.93,80.83,82.38,83.66,84.73,85.64,86.45,87.02,87.51,87.81,88.04])
grod_s30=np.array([-70.41,-69.79,-69.27,-68.78,-68.36,-68.22,-68.51,-69.34,-70.69,-72.56,-74.74,-76.96,-78.98,-80.64,-81.92,-82.87,-83.56,-84.07,-84.4,-84.6,-84.68,-84.65,-84.52,-84.27,-83.87,-83.34,-82.62,-81.7,-80.58,-79.27,-77.82,-76.31,-74.85,-73.45,-72.23,-71.21,-70.41])

longit=360-np.array([0,10,20,30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210,220,230,240,250,260,270,280,290,300,310,320,330,340,350,360])



    # In[2]

sector1_lat=np.array([])
sector1_lon=np.array([])
sector1_loc=np.array([])
sector2_lat=np.array([])
sector2_lon=np.array([])
sector2_loc=np.array([])
sector3_lat=np.array([])
sector3_lon=np.array([])
sector3_loc=np.array([])
sector4_lat=np.array([])
sector4_lon=np.array([])
sector4_loc=np.array([])
sector5_lat=np.array([])
sector5_lon=np.array([])
sector5_loc=np.array([])

sector6_lat=np.array([])
sector6_lon=np.array([])
sector6_loc=np.array([])
sector7_lat=np.array([])
sector7_lon=np.array([])
sector7_loc=np.array([])



sector1_lat150=np.array([])
sector1_lon150=np.array([])
sector1_loc150=np.array([])
sector2_lat150=np.array([])
sector2_lon150=np.array([])
sector2_loc150=np.array([])
sector3_lat150=np.array([])
sector3_lon150=np.array([])
sector3_loc150=np.array([])
sector4_lat150=np.array([])
sector4_lon150=np.array([])
sector4_loc150=np.array([])
sector5_lat150=np.array([])
sector5_lon150=np.array([])
sector5_loc150=np.array([])

sector6_lat150=np.array([])
sector6_lon150=np.array([])
sector6_loc150=np.array([])
sector7_lat150=np.array([])
sector7_lon150=np.array([])
sector7_loc150=np.array([])


# In[5]

#  code that attempts to find the point of infinity in each dataset

cml=100

ofl_lat=np.array([])
ofl_lon=np.array([])

s_ssl=str(int(cml/10+1)*10)

# print(s_ssl)


v = open('fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt')
north_plotting = np.zeros([28,36,4])

header = v.readline()

null1a='0'
null1b='0'
null2a='0'
null2b='0'

n1=np.array([0.])
n1b=np.array([.0])
n2=np.array([0.])
n2b=np.array([.0])

# for z in range(28):
#     for y in range(36):
        
#         line = v.readline()

#         null2a=line[21:29]
#         null2b=line[30:38]
#         if any(chr.isdigit() for chr in null2a): #null2 is real
#             if any(chr.isdigit() for chr in null1a):
#                 donothing=1
#             else: # but null 1 is not real
#                 n2[0]=null2.strip()
#                 n2b[0]=null2b.strip()
#                 ofl_lat=np.append(ofl_lat,n2[0])
#                 ofl_lon=np.append(ofl_lon,n2b[0])
                
#         else: #null2 is not real
#             if any(chr.isdigit() for chr in null1a):# but null 1 is real
#                 n1[0]=null1a.strip()
#                 n1b[0]=null1b.strip()
#                 ofl_lat=np.append(ofl_lat,n1[0])
#                 ofl_lon=np.append(ofl_lon,n1b[0])
 
#         null1a=null2a #null1 is the value in the last position checked
#         null1b=null2b
        
# print(ofl_lat)
# print(ofl_lon)


        
        
# In[4]
        
for cml in range(50,329,10):
    # cml=jc
    # # print(cml)
    # cml=float(80)
    # cml=float(140)
    # # cml=float(180)
    # # cml=float(210)
    # # cml=float(270)
    # cml=float(100)
    js=80
    
    subsolarlat=0
    subsolarlongitude=cml#240.3#+180
    
    
    s_ssl=str(int(subsolarlongitude/10+1)*10)
    
    # print(s_ssl)
    
    
    v = open('fieldline_tracing/fieldline_tracing_mapping_contours_kk2009ext_jrm09int_north_sslong'+s_ssl+'.txt')
    north_plotting = np.zeros([28,36,4])
    
    header = v.readline()
    
    for z in range(28):
        for y in range(36):
            line = v.readline()
            north_plotting[z,y,0]=line[0:12]
            north_plotting[z,y,1]=line[14:20]
            null1=line[21:29]
            if any(chr.isdigit() for chr in null1):
                north_plotting[z,y,2]=null1
            else:
                north_plotting[z,y,2]=np.nan
            null2=line[30:38]
            if any(chr.isdigit() for chr in null2):
                north_plotting[z,y,3]=null2
            else:
                north_plotting[z,y,3]=np.nan
                           
            null2a=line[21:29]
            null2b=line[30:38]
            if any(chr.isdigit() for chr in null2a): #null2 is real
                if any(chr.isdigit() for chr in null1a):
                    donothing=1
                else: # but null 1 is not real
                    n2[0]=null2a.strip()
                    n2b[0]=null2b.strip()
                    # print('****',n2,'****',n2b)
                    if n2[0] < 85:
                        ofl_lat=np.append(ofl_lat,n2[0])
                        ofl_lon=np.append(ofl_lon,n2b[0])
                    
            else: #null2 is not real
                if any(chr.isdigit() for chr in null1a):# but null 1 is real
                    n1[0]=null1a.strip()
                    n1b[0]=null1b.strip()
                    ofl_lat=np.append(ofl_lat,n1[0])
                    ofl_lon=np.append(ofl_lon,n1b[0])
     
            null1a=null2a #null1 is the value in the last position checked
            null1b=null2b

    # calculate the longitude range in this way,              from cml
    
    # for i in range(1,6):
    #     # print(i,(75+i*30-180+360+cml) % 360,(75+(i+1)*30-180+360+cml) % 360)
    #     lon_min=(75+i*30-180+360+cml) % 360
    #     lon_max=(75+(i+1)*30-180+360+cml) % 360
    #     np_15_lon=north_plotting[0,:,3]
    #     np_15_lat=north_plotting[0,:,2]
    #     np_15_loc=north_plotting[0,:,1]
        
    #     np_15_loc=np_15_loc[(np_15_lon < lon_max) & (np_15_lon < lon_min)]
    #     np_15_lat=np_15_lat[(np_15_lon < lon_max) & (np_15_lon < lon_min)]
    #     np_15_lon=np_15_lon[(np_15_lon < lon_max) & (np_15_lon < lon_min)]
    


    ############################################
    ######plotting uncomment
    ############################################
    # crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
    # ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))
    
    # ax.set_global()
    # ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))
    
    
    # gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
    # gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
    # gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)
    # # gl.xformatter = LONGITUDE_FORMATTER
    # # gl.yformatter = LATITUDE_FORMATTER
    
    
    # plt.plot(longit,              grod_n30,color='blue',              linewidth=1,transform=ccrs.PlateCarree())
    # plt.plot(longit,              grod_s30,color='blue',              linewidth=1,transform=ccrs.PlateCarree())
    # plt.plot(longit,              grod_nio,color='darkblue',              linestyle='--',linewidth=1,transform=ccrs.PlateCarree())
    # plt.plot(longit,              grod_sio,color='darkblue',              linestyle='--',linewidth=1,transform=ccrs.PlateCarree())
    # plt.plot([360-185],              [75],color='darkblue',              marker='2',linewidth=1,transform=ccrs.PlateCarree())
    # plt.plot([360-cml,360-cml],              [85,-85],color='darkblue',              linestyle='dotted',linewidth=1,transform=ccrs.PlateCarree())
      
     
    
    # five sectors 
    for i in range(1,6):
        # print(i,(i*30-180+75) ,((i+1)*30-180+75) )
        # print(i,(cml+i*30-180+75) ,(cml+(i+1)*30-180+75) )
        lon_min=360-(((i+1)*30-180+75+cml) % 360)
        lon_max=360-((i*30-180+75+cml) % 360)
        np_15_lon2=north_plotting[0,:,3]
        np_15_lat2=north_plotting[0,:,2]
        np_15_loc2=north_plotting[0,:,1]

        np_150_lon2=north_plotting[27,:,3]
        np_150_lat2=north_plotting[27,:,2]
        np_150_loc2=north_plotting[27,:,1]
    
        # print(np_15_loc.shape)
    
        np_15_loc=np_15_loc2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
        np_15_lat=np_15_lat2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
        np_15_lon=np_15_lon2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
 
    
        np_150_loc=np_150_loc2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]
        np_150_lat=np_150_lat2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]
        np_150_lon=np_150_lon2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]

        # print(np_15_loc.shape)
    
        if i == 1:
            if np_15_lat.size > 0:
                sector1_lat=np.append(sector1_lat,np_15_lat)
                sector1_lon=np.append(sector1_lon,np_15_lon)
                sector1_loc=np.append(sector1_loc,np_15_loc)
            if np_150_lat.size > 0:
                sector1_lat150=np.append(sector1_lat150,np_150_lat)
                sector1_lon150=np.append(sector1_lon150,np_150_lon)
                sector1_loc150=np.append(sector1_loc150,np_150_loc)
        if i == 2:
            if np_15_lat.size > 0:
                sector2_lat=np.append(sector2_lat,np_15_lat)
                sector2_lon=np.append(sector2_lon,np_15_lon)
                sector2_loc=np.append(sector2_loc,np_15_loc)
            if np_150_lat.size > 0:
                 sector2_lat150=np.append(sector2_lat150,np_150_lat)
                 sector2_lon150=np.append(sector2_lon150,np_150_lon)
                 sector2_loc150=np.append(sector2_loc150,np_150_loc)
        if i == 3:
             if np_15_lat.size > 0:
                 sector3_lat=np.append(sector3_lat,np_15_lat)
                 sector3_lon=np.append(sector3_lon,np_15_lon)
                 sector3_loc=np.append(sector3_loc,np_15_loc)
             if np_150_lat.size > 0:
                sector3_lat150=np.append(sector3_lat150,np_150_lat)
                sector3_lon150=np.append(sector3_lon150,np_150_lon)
                sector3_loc150=np.append(sector3_loc150,np_150_loc)
        if i == 4:
             if np_15_lat.size > 0:
                 sector4_lat=np.append(sector4_lat,np_15_lat)
                 sector4_lon=np.append(sector4_lon,np_15_lon)
                 sector4_loc=np.append(sector4_loc,np_15_loc)
             if np_150_lat.size > 0:
                 sector4_lat150=np.append(sector4_lat150,np_150_lat)
                 sector4_lon150=np.append(sector4_lon150,np_150_lon)
                 sector4_loc150=np.append(sector4_loc150,np_150_loc)
        if i == 5:
              if np_15_lat.size > 0:
                  sector5_lat=np.append(sector5_lat,np_15_lat)
                  sector5_lon=np.append(sector5_lon,np_15_lon)
                  sector5_loc=np.append(sector5_loc,np_15_loc)
              if np_150_lat.size > 0:
                  sector5_lat150=np.append(sector5_lat150,np_150_lat)
                  sector5_lon150=np.append(sector5_lon150,np_150_lon)
                  sector5_loc150=np.append(sector5_loc150,np_150_loc)
             

# seven sectors
    for i in range(1,8):
        # print(i,(i*30-180+75) ,((i+1)*30-180+75) )
        # print(i,(cml+i*30-180+75) ,(cml+(i+1)*30-180+75) )
        lon_min=360-(((i+1)*20-180+90+cml) % 360)
        lon_max=360-((i*20-180+90+cml) % 360)
        np_15_lon2=north_plotting[0,:,3]
        np_15_lat2=north_plotting[0,:,2]
        np_15_loc2=north_plotting[0,:,1]
    
    
        np_150_lon2=north_plotting[27,:,3]
        np_150_lat2=north_plotting[27,:,2]
        np_150_loc2=north_plotting[27,:,1]
    
    
        # print(np_15_loc.shape)
    
        np_15_loc=np_15_loc2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
        np_15_lat=np_15_lat2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
        np_15_lon=np_15_lon2[(np_15_lon2 < lon_max) & (np_15_lon2 > lon_min) & (np.abs(np_15_lat2) < 86)]
 
        np_150_loc=np_150_loc2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]
        np_150_lat=np_150_lat2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]
        np_150_lon=np_150_lon2[(np_150_lon2 < lon_max) & (np_150_lon2 > lon_min) & (np.abs(np_150_lat2) < 86)]

    
        # print(np_15_loc.shape)
    
        if i == 1:
            if np_15_lat.size > 0:
                sector1_lat=np.append(sector1_lat,np_15_lat)
                sector1_lon=np.append(sector1_lon,np_15_lon)
                sector1_loc=np.append(sector1_loc,np_15_loc)
            if np_150_lat.size > 0:
                sector1_lat150=np.append(sector1_lat150,np_150_lat)
                sector1_lon150=np.append(sector1_lon150,np_150_lon)
                sector1_loc150=np.append(sector1_loc150,np_150_loc)
        if i == 2:
            if np_15_lat.size > 0:
                sector2_lat=np.append(sector2_lat,np_15_lat)
                sector2_lon=np.append(sector2_lon,np_15_lon)
                sector2_loc=np.append(sector2_loc,np_15_loc)
            if np_150_lat.size > 0:
                 sector2_lat150=np.append(sector2_lat150,np_150_lat)
                 sector2_lon150=np.append(sector2_lon150,np_150_lon)
                 sector2_loc150=np.append(sector2_loc150,np_150_loc)
        if i == 3:
             if np_15_lat.size > 0:
                 sector3_lat=np.append(sector3_lat,np_15_lat)
                 sector3_lon=np.append(sector3_lon,np_15_lon)
                 sector3_loc=np.append(sector3_loc,np_15_loc)
             if np_150_lat.size > 0:
                sector3_lat150=np.append(sector3_lat150,np_150_lat)
                sector3_lon150=np.append(sector3_lon150,np_150_lon)
                sector3_loc150=np.append(sector3_loc150,np_150_loc)
        if i == 4:
             if np_15_lat.size > 0:
                 sector4_lat=np.append(sector4_lat,np_15_lat)
                 sector4_lon=np.append(sector4_lon,np_15_lon)
                 sector4_loc=np.append(sector4_loc,np_15_loc)
             if np_150_lat.size > 0:
                 sector4_lat150=np.append(sector4_lat150,np_150_lat)
                 sector4_lon150=np.append(sector4_lon150,np_150_lon)
                 sector4_loc150=np.append(sector4_loc150,np_150_loc)
        if i == 5:
              if np_15_lat.size > 0:
                  sector5_lat=np.append(sector5_lat,np_15_lat)
                  sector5_lon=np.append(sector5_lon,np_15_lon)
                  sector5_loc=np.append(sector5_loc,np_15_loc)
              if np_150_lat.size > 0:
                  sector5_lat150=np.append(sector5_lat150,np_150_lat)
                  sector5_lon150=np.append(sector5_lon150,np_150_lon)
                  sector5_loc150=np.append(sector5_loc150,np_150_loc)
        if i == 6:
               if np_15_lat.size > 0:
                   sector6_lat=np.append(sector6_lat,np_15_lat)
                   sector6_lon=np.append(sector6_lon,np_15_lon)
                   sector6_loc=np.append(sector6_loc,np_15_loc)
               if np_150_lat.size > 0:
                 sector6_lat150=np.append(sector6_lat150,np_150_lat)
                 sector6_lon150=np.append(sector6_lon150,np_150_lon)
                 sector6_loc150=np.append(sector6_loc150,np_150_loc)
        if i == 7:
              if np_15_lat.size > 0:
                  sector7_lat=np.append(sector7_lat,np_15_lat)
                  sector7_lon=np.append(sector7_lon,np_15_lon)
                  sector7_loc=np.append(sector7_loc,np_15_loc)
              if np_150_lat.size > 0:
                 sector7_lat150=np.append(sector7_lat150,np_150_lat)
                 sector7_lon150=np.append(sector7_lon150,np_150_lon)
                 sector7_loc150=np.append(sector7_loc150,np_150_loc)
             



    # In[4]

cml=180
subsolarlongitude=120
crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector1_lon,sector1_lat,transform=ccrs.PlateCarree(),c=sector1_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
plt.scatter(sector1_lon150,sector1_lat150,transform=ccrs.PlateCarree(),c=sector1_loc150,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.show()


subsolarlongitude=150

crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector2_lon,sector2_lat,transform=ccrs.PlateCarree(),c=sector2_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
plt.scatter(sector2_lon150,sector2_lat150,transform=ccrs.PlateCarree(),c=sector2_loc150,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.show()
subsolarlongitude=180

crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector3_lon,sector3_lat,transform=ccrs.PlateCarree(),c=sector3_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
plt.scatter(sector3_lon150,sector3_lat150,transform=ccrs.PlateCarree(),c=sector3_loc150,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.show()

subsolarlongitude=210

crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector4_lon,sector4_lat,transform=ccrs.PlateCarree(),c=sector4_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
plt.scatter(sector4_lon150,sector4_lat150,transform=ccrs.PlateCarree(),c=sector4_loc150,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.show()

subsolarlongitude=240

crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              js))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector5_lon,sector5_lat,transform=ccrs.PlateCarree(),c=sector5_loc,vmin=0,vmax=24,cmap='hsv',marker='.')
plt.scatter(sector5_lon150,sector5_lat150,transform=ccrs.PlateCarree(),c=sector5_loc150,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.show()
# usage

# In[7]

from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse


crs = ccrs.RotatedPole(globe=ccrs.Globe(flattening=(0.0)))
ax = plt.axes(projection=ccrs.Orthographic(360-cml,              80))

ax.set_global()
ax.add_feature(jup_night(subsolarlat,360-subsolarlongitude,              alpha=0.2))

gl = ax.gridlines(crs=ccrs.PlateCarree(),              linewidth=1,              color='black',              alpha=0.5,              linestyle='dotted',              draw_labels=True)
gl.xlocator = mticker.FixedLocator(np.arange(18)*20-180)
gl.ylocator = mticker.FixedLocator(np.arange(9)*20-90)


plt.scatter(sector5_lon,sector5_lat,transform=ccrs.PlateCarree(),c=sector5_loc,vmin=0,vmax=24,cmap='hsv',marker='.')

plt.scatter(ofl_lon,ofl_lat,marker='.',transform=ccrs.PlateCarree(),linewidth=0.01)


print(ofl_lat.size,ofl_lon.size)


n_elements=ofl_lat.size
a_points = np.zeros([n_elements,2])

a_points[:,0]=ofl_lat
a_points[:,1]=ofl_lon


print(a_points[:,0])

ell = EllipseModel()
ell.estimate(a_points)


xc,              yc,              a,              b,              theta = ell.params

print("center = ",               (xc,              yc))
print("angle of rotation = ",               theta)
print("axes = ",              (a,b))


fig,              axs = plt.subplots(2,              1,              sharex=True,              sharey=True)
axs[0].scatter(a_points[:,0],a_points[:,1])

axs[1].scatter(a_points[:,0],              a_points[:,1])
axs[1].scatter(xc,              yc,              color='red',              s=100)
axs[1].set_xlim(55,              100)
axs[1].set_ylim(120,              200)

ell_patch = Ellipse((xc,              yc),              2*a,              2*b,              theta*180/np.pi,              edgecolor='red',              facecolor='none')

axs[1].add_patch(ell_patch)
plt.show()


# In[8]

sector1_loc[sector1_loc<1]=sector1_loc[sector1_loc<1]+24

plt.plot(sector1_lon,sector1_loc,marker='.',linestyle="None")
plt.plot(sector2_lon,sector2_loc,marker='.',linestyle="None")
plt.plot(sector3_lon,sector3_loc,marker='.',linestyle="None")
plt.plot(sector4_lon,sector4_loc,marker='.',linestyle="None")
plt.plot(sector5_lon,sector5_loc,marker='.',linestyle="None")

plt.plot(sector6_lon,sector6_loc,marker='.',linestyle="None")
plt.plot(sector7_lon,sector7_loc,marker='.',linestyle="None")


x=np.linspace(90,230,140)
poly1 = np.polyfit(sector1_lon,sector1_loc,              deg=5)
poly2 = np.polyfit(sector2_lon,sector2_loc,              deg=5)
poly3 = np.polyfit(sector3_lon,sector3_loc,              deg=5)
poly4 = np.polyfit(sector4_lon,sector4_loc,              deg=5)
poly5 = np.polyfit(sector5_lon,sector5_loc,              deg=5)

poly6 = np.polyfit(sector6_lon,sector6_loc,              deg=5)
poly7 = np.polyfit(sector7_lon,sector7_loc,              deg=5)


plt.plot(x,np.polyval(poly1,              x),              label='fit')
plt.plot(x,np.polyval(poly2,              x),              label='fit')
plt.plot(x,np.polyval(poly3,              x),              label='fit')
plt.plot(x,np.polyval(poly4,              x),              label='fit')
plt.plot(x,np.polyval(poly5,              x),              label='fit')
plt.plot(x,np.polyval(poly6,              x),              label='fit')
plt.plot(x,np.polyval(poly7,              x),              label='fit')
plt.show()

# In[9]

sector1_loc[sector1_loc<1]=sector1_loc[sector1_loc<1]+24

plt.plot(sector1_lon150,sector1_loc150,marker='.',linestyle="None")
plt.plot(sector2_lon150,sector2_loc150,marker='.',linestyle="None")
plt.plot(sector3_lon150,sector3_loc150,marker='.',linestyle="None")
plt.plot(sector4_lon150,sector4_loc150,marker='.',linestyle="None")
plt.plot(sector5_lon150,sector5_loc150,marker='.',linestyle="None")

plt.plot(sector6_lon150,sector6_loc150,marker='.',linestyle="None")
plt.plot(sector7_lon150,sector7_loc150,marker='.',linestyle="None")


x=np.linspace(120,230,110)
poly1 = np.polyfit(sector1_lon150,sector1_loc150,              deg=5)
poly2 = np.polyfit(sector2_lon150,sector2_loc150,              deg=5)
poly3 = np.polyfit(sector3_lon150,sector3_loc150,              deg=5)
poly4 = np.polyfit(sector4_lon150,sector4_loc150,              deg=5)
poly5 = np.polyfit(sector5_lon150,sector5_loc150,              deg=5)

poly6 = np.polyfit(sector6_lon150,sector6_loc150,              deg=5)
poly7 = np.polyfit(sector7_lon150,sector7_loc150,              deg=5)


# plt.plot(x,np.polyval(poly1,              x),              label='fit')
# plt.plot(x,np.polyval(poly2,              x),              label='fit')
# plt.plot(x,np.polyval(poly3,              x),              label='fit')
# plt.plot(x,np.polyval(poly4,              x),              label='fit')
# plt.plot(x,np.polyval(poly5,              x),              label='fit')
# plt.plot(x,np.polyval(poly6,              x),              label='fit')
# plt.plot(x,np.polyval(poly7,              x),              label='fit')
# plt.yrange=[0,24]
plt.show()





# jdiam,jcml,jsel,jd,ja =   get_values('05jun99')   
    
# print(jsel)

